The SMTK contains two modules for the characterisation of ground motion:
1) smtk.response_spectrum
This module contains methods for calculation of the response of a set of single degree-of-freedom (SDOF) oscillators using an input time series. Two methods are currently supported:
i) Newmark-Beta
ii) Nigam & Jennings (1969) {Preferred}
The module also includes functions for plotting the response spectra and time series
2) smtk.intensity_measures
This module contains a set of functions for deriving different intensity measures from a strong motion record
i) get_peak_measures(...) - returns PGA, PGV and PGD
ii) get_response_spectrum(...) - returns the response spectrum
iii) get_response_spectrum_pair(...) - returns a response spectrum pair
iv) geometric_mean_spectrum(...) - returns the geometric mean of a pair of records
v) arithmetic_mean_spectrum(...) - returns the arithmetic mean of a pair of records
vi) geometric_mean_spectrum(...) - returns the envelope spectrum of a pair of records
vii) larger_pga(...) - Returns the spectrum with the larger PGA
viii) rotate_horizontal(...) - rotates a record pair through angle theta
ix) gmrotdpp(...) - Returns the rotationally-dependent geometric fractile (pp) of a pair of records
x) gmrotipp(...) - Returns the rotationally-independent geometric fractile (pp) of a pair of records
In [ ]:
%load_ext autoreload
%autoreload 2
import warnings; warnings.filterwarnings("ignore")
In [ ]:
# Import modules
%matplotlib inline
import numpy as np # Numerical Python package
import matplotlib.pyplot as plt # Python plotting package
# Import
import smtk.response_spectrum as rsp # Response Spectra tools
import smtk.intensity_measures as ims # Intensity Measure Tools
periods = np.array([0.01, 0.02, 0.03, 0.04, 0.05, 0.075, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19,
0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.44, 0.46, 0.48, 0.5,
0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.5, 6.0,
6.5, 7.0,7.5, 8.0, 8.5, 9.0, 9.5, 10.0], dtype=float)
number_periods = len(periods)
In [ ]:
# Load record pair from files
x_record = np.genfromtxt("data/sm_record_x.txt")
y_record = np.genfromtxt("data/sm_record_y.txt")
x_time_step = 0.002 # Record sampled at 0.002 s
y_time_step = 0.002
In [ ]:
# Create an instance of the Newmark-Beta class
nigam_jennings = rsp.NigamJennings(x_record, x_time_step, periods, damping=0.05, units="cm/s/s")
sax, time_series, acc, vel, dis = nigam_jennings.evaluate()
# Plot Response Spectrum
rsp.plot_response_spectra(sax, axis_type="semilogx", filename="images/response_nigam_jennings.pdf",filetype="pdf")
In [ ]:
rsp.plot_time_series(time_series["Acceleration"],
x_time_step,
time_series["Velocity"],
time_series["Displacement"])
In [ ]:
pga_x, pgv_x, pgd_x, _, _ = ims.get_peak_measures(0.002, x_record, True, True)
print "PGA = %10.4f cm/s/s, PGV = %10.4f cm/s, PGD = %10.4f cm" % (pga_x, pgv_x, pgd_x)
pga_y, pgv_y, pgd_y, _, _ = ims.get_peak_measures(0.002, y_record, True, True)
print "PGA = %10.4f cm/s/s, PGV = %10.4f cm/s, PGD = %10.4f cm" % (pga_y, pgv_y, pgd_y)
In [ ]:
print "Bracketed Duration (> 5 cm/s/s) = %9.3f s" % ims.get_bracketed_duration(x_record, x_time_step, 5.0)
print "Uniform Duration (> 5 cm/s/s) = %9.3f s" % ims.get_uniform_duration(x_record, x_time_step, 5.0)
print "Significant Duration (5 - 95 Arias ) = %9.3f s" % ims.get_significant_duration(x_record, x_time_step, 0.05, 0.95)
In [ ]:
print "Arias Intensity = %12.4f cm-s" % ims.get_arias_intensity(x_record, x_time_step)
print "Arias Intensity (5 - 95) = %12.4f cm-s" % ims.get_arias_intensity(x_record, x_time_step, 0.05, 0.95)
print "CAV = %12.4f cm-s" % ims.get_cav(x_record, x_time_step)
print "CAV5 = %12.4f cm-s" % ims.get_cav(x_record, x_time_step, threshold=5.0)
print "Arms = %12.4f cm-s" % ims.get_arms(x_record, x_time_step)
In [ ]:
# Get response spectrum
sax = ims.get_response_spectrum(x_record, x_time_step, periods)[0]
print "Velocity Spectrum Intensity (cm/s/s) = %12.5f" % ims.get_response_spectrum_intensity(sax)
print "Acceleration Spectrum Intensity (cm-s) = %12.5f" % ims.get_acceleration_spectrum_intensity(sax)
In [ ]:
sax, say = ims.get_response_spectrum_pair(x_record, x_time_step,
y_record, y_time_step,
periods,
damping=0.05,
units="cm/s/s",
method="Nigam-Jennings")
In [ ]:
sa_gm = ims.geometric_mean_spectrum(sax, say)
rsp.plot_response_spectra(sa_gm, "semilogx", filename="images/geometric_mean_spectrum.pdf", filetype="pdf")
In [ ]:
sa_env = ims.envelope_spectrum(sax, say)
rsp.plot_response_spectra(sa_env, "semilogx", filename="images/envelope_spectrum.pdf", filetype="pdf")
In [ ]:
gmrotd50 = ims.gmrotdpp(x_record, x_time_step, y_record, y_time_step, periods, percentile=50.0,
damping=0.05, units="cm/s/s")
gmroti50 = ims.gmrotipp(x_record, x_time_step, y_record, y_time_step, periods, percentile=50.0,
damping=0.05, units="cm/s/s")
In [ ]:
# Plot all of the rotational angles!
plt.figure(figsize=(8, 6))
for row in gmrotd50["GeoMeanPerAngle"]:
plt.semilogx(periods, row, "-", color="LightGray")
plt.semilogx(periods, gmrotd50["GMRotDpp"], 'b-', linewidth=2, label="GMRotD50")
plt.semilogx(periods, gmroti50["Pseudo-Acceleration"], 'r-', linewidth=2, label="GMRotI50")
plt.xlabel("Period (s)", fontsize=18)
plt.ylabel("Acceleration (cm/s/s)", fontsize=18)
plt.legend(loc=0)
plt.savefig("images/rotational_spectra.pdf", dpi=300, format="pdf")
In [ ]:
ims.plot_fourier_spectrum(x_record, x_time_step,
filename="images/fourier_spectrum.pdf", filetype="pdf")
In [ ]:
from smtk.smoothing.konno_ohmachi import KonnoOhmachi
# Get the original Fourier spectrum
freq, amplitude = ims.get_fourier_spectrum(x_record, x_time_step)
# Configure Smoothing Parameters
smoothing_config = {"bandwidth": 40, # Size of smoothing window (lower = more smoothing)
"count": 1, # Number of times to apply smoothing (may be more for noisy records)
"normalize": True}
# Apply the Smoothing
smoother = KonnoOhmachi(smoothing_config)
smoothed_spectra = smoother.apply_smoothing(amplitude, freq)
# Compare the Two Spectra
plt.figure(figsize=(7,5))
plt.loglog(freq, amplitude, "k-", lw=1.0,label="Original")
plt.loglog(freq, smoothed_spectra, "r", lw=2.0, label="Smoothed")
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.xlim(0.05, 200)
plt.ylabel("Fourier Amplitude", fontsize=14)
plt.tick_params(labelsize=12)
plt.legend(loc=0, fontsize=14)
plt.grid(True)
plt.savefig("images/SmoothedFourierSpectra.pdf", format="pdf", dpi=300)
In [ ]:
# Load in a three component data set
record_file = "data/record_3component.csv"
record_3comp = np.genfromtxt(record_file, delimiter=",")
time_vector = record_3comp[:, 0]
x_record = record_3comp[:, 1]
y_record = record_3comp[:, 2]
v_record = record_3comp[:, 3]
time_step = 0.002
# Plot the records
fig = plt.figure(figsize=(8,12))
fig.set_tight_layout(True)
ax = plt.subplot(311)
ax.plot(time_vector, x_record)
ax.set_ylim(-80., 80.)
ax.set_xlim(0., 10.5)
ax.grid(True)
ax.set_xlabel("Time (s)", fontsize=14)
ax.set_ylabel("Acceleration (cm/s/s)", fontsize=14)
ax.tick_params(labelsize=12)
ax.set_title("EW", fontsize=16)
ax = plt.subplot(312)
ax.plot(time_vector, y_record)
ax.set_xlim(0., 10.5)
ax.set_ylim(-80., 80.)
ax.grid(True)
ax.set_xlabel("Time (s)", fontsize=14)
ax.set_ylabel("Acceleration (cm/s/s)", fontsize=14)
ax.set_title("NS", fontsize=16)
ax.tick_params(labelsize=12)
ax = plt.subplot(313)
ax.plot(time_vector, v_record)
ax.set_xlim(0., 10.5)
ax.set_ylim(-40., 40.)
ax.grid(True)
ax.set_xlabel("Time (s)", fontsize=14)
ax.set_ylabel("Acceleration (cm/s/s)", fontsize=14)
ax.set_title("Vertical", fontsize=16)
ax.tick_params(labelsize=12)
plt.savefig("images/3component_timeseries.pdf", format="pdf", dpi=300)
In [ ]:
x_freq, x_four = ims.get_fourier_spectrum(x_record, time_step)
y_freq, y_four = ims.get_fourier_spectrum(y_record, time_step)
v_freq, v_four = ims.get_fourier_spectrum(v_record, time_step)
plt.figure(figsize=(7, 5))
plt.loglog(x_freq, x_four, "k-", lw=1.0, label="EW")
plt.loglog(y_freq, y_four, "b-", lw=1.0, label="NS")
plt.loglog(v_freq, v_four, "r-", lw=1.0, label="V")
plt.xlim(0.05, 200.)
plt.tick_params(labelsize=12)
plt.grid(True)
plt.xlabel("Frequency (Hz)", fontsize=16)
plt.ylabel("Fourier Amplitude", fontsize=16)
plt.legend(loc=3, fontsize=16)
plt.savefig("images/3component_fas.pdf", format="pdf", dpi=300)
In [ ]:
# Setup parameters
params = {"Function": "KonnoOhmachi",
"bandwidth": 40.0,
"count": 1.0,
"normalize": True
}
# Returns
# 1. Horizontal to Vertical Spectral Ratio
# 2. Frequency
# 3. Maximum H/V
# 4. Period of Maximum H/V
hvsr, freq, max_hv, t_0 = ims.get_hvsr(x_record, time_step, y_record, time_step, v_record, time_step, params)
In [ ]:
plt.figure(figsize=(7,5))
plt.semilogx(freq, hvsr, 'k-', lw=2.0)
# Show T0
t_0_line = np.array([[t_0, 0.0],
[t_0, 1.1 * max_hv]])
plt.semilogx(1.0 / t_0_line[:, 0], t_0_line[:, 1], "r--", lw=1.5)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("H / V", fontsize=14)
plt.tick_params(labelsize=14)
plt.xlim(0.1, 10.0)
plt.grid(True)
plt.title(r"$T_0 = %.4f s$" % t_0, fontsize=16)
plt.savefig("images/hvsr_example1.pdf", format="pdf", dpi=300)